function [ results_fnames ] = osl_run_first_level_epoched( oat )

% [ results_fnames ] = osl_run_first_level_epoched( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level trialwise GLM
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2011

global OSLDIR;

first_level=oat.first_level;
source_recon_name='source_recon';

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if(~isfield(first_level,'bc_trialwise')),
    first_level.bc_trialwise=0;
end;

% check if user has asked to jailbreak the trialwise data
if(~isfield(first_level,'save_data'));
    first_level.save_data = 0;
end

if(~isfield(first_level, 'trialwise_directory'));
    first_level.trialwise_directory = [];
end

% set the downsample factor.  If doing a tf, read tf_downsample_factor.  If
% doing time domain, read time_downsample_factor
switch first_level.tf_method
    case 'none' % use time_downsample_factor
        if ~isempty(first_level.time_downsample_factor)
            ds_factor = 1./first_level.time_downsample_factor;
        else
            ds_factor = [];
        end
    otherwise  % use tf_downsample_factor
        if ~isempty(first_level.time_downsample_factor)
            warning('first_level.time_downsample_factor is IGNORED with a time-freq analysis of epoched data. Use first_level.tf_downsample_factor instead.');
            
        end;
        
        if ~isempty(first_level.tf_downsample_factor)
            ds_factor = 1./first_level.tf_downsample_factor;
        else
            ds_factor = [];
        end
end % switch first_level.tf_method

if ~isempty(first_level.time_moving_av_win_size)
    warning('first_level.time_moving_av_win_size is NOT used with epoched data');
end;

contrast_list=first_level.contrast;
for c=1:length(contrast_list),
    contrast_list{c}=contrast_list{c}(:);
end;

clear stats;

if(sum(first_level.bc)>0 && first_level.baseline_timespan(1)>first_level.baseline_timespan(2)),
    error('first_level.baseline_timespan(1)>first_level.baseline_timespan(2)');
end;

if(~isfield(oat.first_level,'design_matrix') && ~isfield(oat.first_level,'design_matrix_summary')),
    error('Design matrix is not specified');
end;

for subi_todo=1:length(first_level.sessions_to_do),   
        
    sub=first_level.sessions_to_do(subi_todo);
    
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this session
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D=source_recon_results.BF.data.D;
    if(~strcmp([oat.source_recon.dirname],D.path)),        
        D=path(D,[oat.source_recon.dirname]);
        D.save;
        D=spm_eeg_load([D.path '/' D.fname]);
        osl_save_oat_results(oat,source_recon_results);    
    end;
    
    %% setup things that are common to all sessions
    if(subi_todo==1),
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            disp('Working in sensor space');
            
            modality='MEG';
            chanindmeg = strmatch(modality, D.chantype);
            first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
            
        else,
            
            % setup std space brain
            first_level_results.gridstep=source_recon_results.gridstep;
            
            S=[];
            S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
            S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
            S.current_level=first_level;
            S.lower_level_mni_coord=source_recon_results.mni_coord;
            S.lower_level_gridstep=source_recon_results.gridstep;
            
            [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
            clear S;
            
        end;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup HMM 
        if strcmp(source_recon_results.recon_method,'none'),
            NK=1;
        else
            NK=numel(source_recon_results.BF.inverse.W.MEG); 
        end;
        
        %% set up potential time frequency bases in tf_settings, note that any time or
        % frequency averaging will be done AFTER the time-frequency decomposition
        
        % set the frequency range
        freq_range = first_level.tf_freq_range; % range used for stat testing
        if isempty(freq_range)
            % if not set as a first_level option, default to the range used
            % for source_recon
            freq_range=source_recon_results.source_recon.freq_range;
            
            if isempty(freq_range)
                freq_range=[1 D.fsample/3];
                disp(['WARNING: No frequency range has been specified. Set oat.source_recon.freq_range or oat.first_level.tf_freq_range. Using ' num2str(freq_range) ' Hz.']);                
            end % isempty(freq_range)
        end % if isempty(freq_range)
        
        %source_recon_time_indices = find(D.time>=source_recon_results.source_recon.time_range(1) & D.time<=source_recon_results.source_recon.time_range(2));
        source_recon_time_indices=find(source_recon_results.samples2use);
        source_recon_times = D.time(source_recon_time_indices);
        
        if isempty(first_level.time_range)
            first_level.time_range=source_recon_results.woi;
        end
        
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.tf_multitaper_ncycles   = first_level.tf_multitaper_ncycles;
        tf_settings.raw_times               = source_recon_times;
        tf_settings.ds_factor               = ds_factor;
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper     = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin      = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth= first_level.tf_multitaper_freqsmooth;
        
        out = osl_tf_transform( tf_settings ); % passed without a data matrix, osl_tf_transform just gets container sizes
        
        if isfield(out, 'tf_morlet_basis')
            first_level.tf_morlet_basis = out.tf_morlet_basis;
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup time and frequency windows/indices for using later
                       
        % setup result time dimension        
        first_level_time_indices  = ( out.tf_times >= first_level.time_range(1) ) & ( out.tf_times <= first_level.time_range(2) );
        first_level_results.times = out.tf_times(first_level_time_indices);                
        if first_level.time_average
            first_level_results.times=mean(first_level_results.times);
            tres=1;
        end
        
        % setup result freq dimension        
        first_level_results.frequencies = out.tf_freqs;
        if first_level.freq_average
            first_level_results.frequencies=mean(first_level_results.frequencies);
        end
        
        % setup result spatial dimension
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);        
        if(first_level.space_average)
            Nvoxels_out=1;
            if(isfield(first_level_results,'mni_coord')),
                first_level_results.mni_coord=mean(first_level_results.mni_coord,1);
            end;
        end
        
        clear out
        
        %% create results containers
        
        % get dimensions:
        ntpts                       = length(first_level_results.times);
        nfreqs                      = length(first_level_results.frequencies);
        
        % create cope and stdcope containers:
        if(first_level.doGLM)
            first_level_results.stdcope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            first_level_results.cope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            varcope_by_state=zeros([ntpts,length(first_level.contrast),nfreqs,NK]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            
            first_level_results.stdcope_by_state=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs,NK]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            first_level_results.cope_by_state=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs,NK]);% num_voxels x num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
        end;
        
        first_level_results.pseudo_zstat_var = zeros([Nvoxels_out,1]);
        
    end; % end of if(sub==1)
    
    %% setup the GLM design matrix for this session
    if(isfield(first_level,'design_matrix')),
        
        first_level.design_matrix;
        
        Ntrials   = size(first_level.design_matrix,1);
        
    else,
        
        % setup trial types
        trialtypes=[];
        triallist=[];

        for i=1:length(source_recon_results.source_recon.conditions), % indexes conditions/triggers within subset
                               
            trigname=source_recon_results.source_recon.conditions{i};
            Ntrialspercond=length(D.pickconditions(trigname)); %% number of trials for this condition
            triallist=[triallist , D.pickconditions(trigname)];           
            trialtypes=[trialtypes ; ones(Ntrialspercond,1).*i];
            
        end;
        
        S3=[];
        S3.trialtypes=trialtypes;
        
        if numel(first_level.design_matrix_summary)==length(oat.source_recon.results_fnames) &&(ischar(first_level.design_matrix_summary{sub})),
            
            S3.Xsummary=first_level.design_matrix_summary{sub};
            
            if(isfield(first_level,'trial_rejects'))
                S3.trial_rejects=first_level.trial_rejects;
            else
                warning('oat.first_level.trial_rejects not set');
            end;
        else
            S3.Xsummary=first_level.design_matrix_summary;
        end;
        
        first_level.x=setup_beamformer_designmatrix(S3);
        
        Ntrials   = length(triallist);
        
    end;
        
    %% do everything from now on on one voxel at a time (for the sake of
    %% RAM) in the first level mask space
    
    if first_level.save_data
        % transform-specific parameters
        switch first_level.tf_method
            case 'none'
                trlwise.type = 'time_domain';
            otherwise
                trlwise.type    = first_level.tf_method;
                trlwise.dat     = nan(Ntrials,ntpts,nfreqs); % if not a tf, the last dimenstion will be vestigial
                trlwise.freqs   = first_level_results.frequencies;
                trlwise.times   = first_level_results.times;
        end
        
        % spatial parameters
        trlwise.mask        = read_avw([oat.source_recon.dirname '/' oat.first_level.name '_mask']);
        trlwise.mnicoords   = first_level_results.mni_coord;
        trlwise.gridstep    = first_level_results.gridstep;
    end
    
    %%%%%%%%%%%%%%%
    %% setup channels
    modality='MEG';
    chanindmeg = strmatch(modality, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors 
        first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
        chanind = chanindmeg;
    end; 
    
    %%%%%%%%%%%%%%%
    %% load in sensor data
    first_level_results.D_sensor_data=D;
    first_level_results.chanind=chanind;
    
    sensor_data=D(chanind, source_recon_time_indices, triallist);
    class_prs=zeros(NK, length(source_recon_time_indices), length(triallist));
    
    classchanind=find(strcmp(D.chanlabels,'Class'));
    if(isempty(classchanind)),
        error(['No ''Class'' chanlabel in: ' D.fname]);
    else
        for si=1:NK,        
            class_samples_inds{si} = (D(classchanind,source_recon_time_indices, triallist)==si);

            state_tinds=(class_samples_inds{si}==1);                                                

            total_time=(source_recon_times(end)-source_recon_times(1))*size(state_tinds,3);
            state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
            txt=['State ' num2str(si) ' is active for ' num2str(state_time) 'secs'];

            if(state_time<10),
                warning([txt '. It will be excluded from the GLM analysis.']);
            else
                disp(txt);
            end;
            
            classpr_chanind=find(strcmp(D.chanlabels,['ClassPr' num2str(si)]));
            class_prs(si,:,:) = D(classpr_chanind,source_recon_time_indices, triallist);

        end;
    end;
    
    %%%%%%%%%%%%%%%%
    %% do everything from now on one voxel at a time (for the sake of
    %% RAM) in the first level mask space                
    if(first_level.doGLM)
        disp('First level COPEs outputted will have dimension Nvoxels x Ntpts x Ncontrasts x Nfreqs:');
        tmp=[size(first_level_results.cope) 1 1];
        disp([num2str(tmp(1:4))]);
    end;
    
    %%%%%%%%%%%%%%%%
    %% now do space averaging 
    if(first_level.space_average),
        error('Pre-first_level space averaging not implemented');
        
        dat=zeros(Ntrials,length(source_recon_time_indices));
        wnorms=zeros(Ntrials,length(source_recon_time_indices));            
        num=0;
        
        S2=[];
        S2.source_recon_results=source_recon_results;
        S2.chanind=chanind;
        S2.class_samples_inds=class_samples_inds;
        S2.class_prs=class_prs;
        S2.use_class_probs=0;

        S2.sensor_data=sensor_data;
        for indind=1:length(first_level_results.mask_indices_in_source_recon), % indexes brain space
            
            weights_usable=0;
            
            S2.voxind=first_level_results.mask_indices_in_source_recon(indind);           

            [ datnew wnormsnew ] = get_voxel_recon_timecourse( S2 );
                                   
            if sum(isnan(squash(dat)))==0,
                num=num+1;
                dat=dat+datnew;
                wnorms=wnorms+wnormsnew;
            end;
            
        end;
        
        dat=dat/num;
        wnorms=wnorms/num;
    end; 
    
    %%%%%%%%%%%%%%%%
    %% setup design matrices that can be reused over voxels
    ntpts=size(state_tinds,2);
    x_stored=cell(NK,ntpts);
    pinvxtx_stored=cell(NK,ntpts);
    pinvx_stored=cell(NK,ntpts);
    contrast_list_mat=cell2mat(contrast_list);
        
    for si=1:NK,
        state_tinds=(class_samples_inds{si}==1);                                                

        total_time=(source_recon_times(end)-source_recon_times(1))*size(state_tinds,3);
        state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));

        if(state_time>10),
            xc=[];
            contrast_list_matc=[];
            for t= 1 : ntpts, % indexes time points within trials

                x=first_level.x;
                x=x(permute(state_tinds(1,t,:),[3 1 2]),:);                        

                if(oat.first_level.do_glm_demean)
                    for pp=1:size(x,2),
                        x(:,pp)=x(:,pp)-mean(x(:,pp));
                    end;
                end;

                x_stored{si,t}=x;
              
                pinvxtx_stored{si,t}=pinv(x'*x);
                pinvx_stored{si,t}=pinv(x);
                
                xc=blkdiag(xc,x);
                 
                contrast_list_matc=blkdiag(contrast_list_matc,contrast_list_mat);
            end;
            contrast_list_matc=sparse(contrast_list_matc);

            xc_stored{si}=sparse(xc);
            pinvxctxc_stored{si}=sparse(pinv(xc'*xc));
            pinvxc_stored{si}=sparse(pinv(xc));
        end;
    end;
    
    
    %%%%%%%%%%%%%%%%
    %% loop over voxels
    ft_progress('init', 'etf');
    S2=[];
    S2.source_recon_results=source_recon_results;
    S2.chanind=chanind;
    S2.class_samples_inds=class_samples_inds;
    S2.sensor_data=sensor_data;
    S2.class_prs=class_prs;
    S2.use_class_probs=0;
        
    clear sensor_data;
    
    for indind=1:Nvoxels_out, % indexes brain space
        
        ft_progress(indind/Nvoxels_out);
               
        if(~first_level.space_average),
            S2.voxind=first_level_results.mask_indices_in_source_recon(indind);

            [ dat wnorms wnorm wnorms_nai wnorm_nai ] = get_voxel_recon_timecourse( S2 );
              
            dat=dat/sqrt(mean(squash(wnorms_nai)));
            %dat=dat./sqrt(wnorms_nai);
            
            first_level_results.pseudo_zstat_var(indind)=var(squash(dat));
         end;
        
        if sum(isnan(squash(dat)))~=0,
            
            if(first_level.doGLM)
                first_level_results.cope(indind,:,:,:)=nan;
                first_level_results.stdcope(indind,:,:,:)=nan;         
                first_level_results.cope_by_state(indind,:,:,:,:)=nan;
                first_level_results.stdcope_by_state(indind,:,:,:,:)=nan; 
            end;
            warning('NANs in recon data');
            
        else
            
            %% Transform the data - either downsampling in time domain, or a TF transform
            
            tf_settings = [];
            tf_settings.tf_method               = first_level.tf_method;
            tf_settings.tf_logtransform         = first_level.tf_logtransform;
            tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
            tf_settings.tf_freq_range           = freq_range; % set above
            tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
            tf_settings.tf_multitaper_ncycles   = first_level.tf_multitaper_ncycles;
            tf_settings.raw_times               = source_recon_times;
            tf_settings.ds_factor               = ds_factor;
            tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
            tf_settings.tf_multitaper_taper     = first_level.tf_multitaper_taper;
            tf_settings.tf_multitaper_twin      = first_level.tf_multitaper_twin;
            tf_settings.tf_multitaper_freqsmooth= first_level.tf_multitaper_freqsmooth;

            if isfield(first_level,'tf_morlet_basis')
                tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
            end
                        
            out = osl_tf_transform( tf_settings , dat ); % do the transformation
            dat = out.dattf;
            
            if first_level.save_data
                trlwise.dat(indind,:,:,:) = dat;
            end % if first_level.save_data
            
            if(first_level.bc_trialwise),
                baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);
                
                if(sum(baseline_time_indices)>0)
                    % normalise by noise in baseline window:
                    for tri=1:size(dat,1), % indexes trials
                        for f = 1 : size(dat,3),
                            dat(tri,:,f)=(dat(tri,:,f)-mean(dat(tri,baseline_time_indices,f)));
                        end;
                    end;
                else
                    warning('Could not baseline correct as there are no timepoints in the specified baseline period');
                end;
            end;
            
            %% now do windowing over time
            dat=dat(:,first_level_time_indices,:);
           
            %% now do averaging over time and frequency
            if(first_level.time_average)
                dat=mean(dat,2);                
            end;
            
            if(first_level.freq_average)
                dat=mean(dat,3);
            end;
            
            %% Do GLM
            if first_level.doGLM,
                if(size(dat,1)>1)                                          

                    for si=1:NK,
                        state_tinds=(class_samples_inds{si}==1);                                                
                              
                        total_time=(source_recon_times(end)-source_recon_times(1))*size(state_tinds,3);
                        state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
                                              
                        if(state_time<10),
                            
                            first_level_results.cope_by_state(indind,:,:,:,si)=0;
                            first_level_results.stdcope_by_state(indind,:,:,:,si)=1e32;
                                    
                        else
                            
                            if(1)

                                for f = 1 : size(dat,3),
                                    dattfc=[];
                                    for t= 1 : size(dat,2), % indexes time points within trials
                                        
                                        dattf=dat(permute(state_tinds(1,t,:),[3 1 2]),t,f);
                                        
                                        if(oat.first_level.do_glm_demean)
                                            dattf=dattf-mean(dattf);                                                        
                                        end;
                                        
                                        dattfc=cat(1,dattfc,dattf);
                                        
                                    end;
                                    
                                    x=xc_stored{si};
                                    pinvxtx=pinvxctxc_stored{si};
                                    pinvx=pinvxc_stored{si};
                                
                                    [copeout, varcopeout, coapeout, dof,pe]=glm_fast_for_meg(dattfc,x,pinvxtx,pinvx,contrast_list_matc,0);
                                    
                                    if(strcmp(oat.first_level.cope_type,'coape')),
                                        first_level_results.cope_by_state(indind,:,:,f,si)=reshape(coapeout,size(contrast_list_mat,2),size(dat,2))';
                                    elseif(strcmp(oat.first_level.cope_type,'cope')),
                                        first_level_results.cope_by_state(indind,:,:,f,si)=reshape(copeout,size(contrast_list_mat,2),size(dat,2))'
                                    elseif(strcmp(oat.first_level.cope_type,'acope')),
                                        first_level_results.cope_by_state(indind,:,:,f,si)=abs(reshape(copeout,size(contrast_list_mat,2),size(dat,2))');
                                        %first_level_results.cope_by_state(indind,:,:,f,si)=abs(reshape(copeout,size(dat,2),size(contrast_list_mat,2)));
                                    else
                                        error('Invalid cope_type');
                                    end;

                                    varcope_by_state(:,:,f,si)=reshape(varcopeout,size(contrast_list_mat,2),size(dat,2))'; 

                                end;                                
                                
                            else
                                
                                error('Not implemented');
                                
                                for t= 1 : size(dat,2), % indexes time points within trials
                                    x=x_stored{si,t};
                                    pinvxtx=pinvxtx_stored{si,t};
                                    pinvx=pinvx_stored{si,t};

                                    for f = 1 : size(dat,3),
                                        dattf=dat(permute(state_tinds(1,t,:),[3 1 2]),t,f);             

                                        if(oat.first_level.do_glm_demean)
                                            dattf=dattf-mean(dattf);                                                        
                                        end;

                                        if(first_level.use_robust_glm),
                                            [b,stats2]=robustfit(x,dattf,'bisquare',4.685,'off'); % fits GLM robustly, Y=XB+E
                                            for c=1:length(contrast_list),

                                                if(strcmp(oat.first_level.cope_type,'coape')),
                                                    first_level_results.cope_by_state(indind,t,c,f,si)=contrast_list{c}'*abs(b);
                                                elseif(strcmp(oat.first_level.cope_type,'cope')),
                                                    first_level_results.cope_by_state(indind,t,c,f,si)=contrast_list{c}'*b;
                                                elseif(strcmp(oat.first_level.cope_type,'acope')),
                                                    first_level_results.cope_by_state(indind,t,c,f,si)=abs(contrast_list{c}'*b);
                                                else
                                                    error('Invalid cope_type');
                                                end;
                                                first_level_results.stdcope_by_state(indind,t,c,f,si)=sqrt(contrast_list{c}'*stats2.covb*contrast_list{c});
                                            end;
                                            pe(t,:)=b; % num_timepoints x  num_pes
                                        else
                                            [copeout, varcopeout, coapeout]=glm_fast_for_meg(dattf,x,pinvxtx,pinvx,contrast_list_mat,0);
                                            if(strcmp(oat.first_level.cope_type,'coape')),
                                                first_level_results.cope_by_state(indind,t,:,f,si)=coapeout;
                                            elseif(strcmp(oat.first_level.cope_type,'cope')),
                                                first_level_results.cope_by_state(indind,t,:,f,si)=copeout;
                                            elseif(strcmp(oat.first_level.cope_type,'acope')),
                                                first_level_results.cope_by_state(indind,t,:,f,si)=abs(copeout);
                                            else
                                                error('Invalid cope_type');
                                            end;

                                            varcope_by_state(t,:,f,si)=varcopeout;
                                            %pe(t,:)=peout; % num_timepoints x  num_pes
                                        end;

                                    end;

                                end;
                            end;
                        end; % si (state index) 
                    end;
                    
                    %%
                    if(NK>1)
                        %% do fixed effects averaging of states
                        for iContrast = 1:size(first_level_results.cope_by_state,3)
                            % freq loop
                            for iFreq = 1:size(first_level_results.cope_by_state,4)
                                % time loop
                                for iTime = 1 : size(first_level_results.cope_by_state,2), % indexes time points within trials
                        
                                    cope_tmp=permute(first_level_results.cope_by_state(indind,iTime,iContrast,iFreq,:),[5 1 2 3 4]);
                                    varcope_tmp=permute(varcope_by_state(iTime,iContrast,iFreq,:),[4 1 2 3 ]);
                                    
                                    sinds=find(varcope_tmp>0);
                                    
                                    z=ones(length(sinds),1);
                                    
                                    cope_tmp=cope_tmp(sinds);
                                    varcope_tmp=varcope_tmp(sinds);
                                    
                                    S=varcope_tmp;
                                    
                                    [gam, ~, covgam] = flame1(cope_tmp,z,S,1);
                                    gam=mean(cope_tmp);
                                    
                                    first_level_results.cope(indind,iTime,iContrast,iFreq)=gam;
                                    first_level_results.stdcope(indind,iTime,iContrast,iFreq)=sqrt(covgam);

                                    % iTime=86;iContrast=2;iFreq=2;figure;plot(squeeze(first_level_results.cope_by_state(indind,:,iContrast,iFreq,:)))
                                   
                                end; % for iTime
                            end % for iFreq = 1:nFreqs
                        end % for iContrast = 1:nContrasts
                         if(1)
                             indind
                             iContrast=[2,4,5,6];iFreq=2;
                             first_level_results.cope(indind,84:88,iContrast,iFreq)
                             first_level_results.stdcope(indind,84:88,iContrast,iFreq)
                             first_level_results.cope(indind,84:88,iContrast,iFreq)./first_level_results.stdcope(indind,84:88,iContrast,iFreq)
                             figure(1);iFreq=2;plot(squeeze(first_level_results.cope(indind,:,iContrast,iFreq)))
                             figure(2);iContrast=2;iFreq=2;plot(squeeze(first_level_results.cope_by_state(indind,:,iContrast,iFreq,:)))
                             figure(3);iContrast=4;iFreq=2;plot(squeeze(first_level_results.cope_by_state(indind,:,iContrast,iFreq,:)))
                             
                             
                             
                             figure(2);iFreq=2;plot(squeeze(first_level_results.stdcope(indind,:,iContrast,iFreq)))
                             figure(4);iContrast=2;iFreq=2;plot(squeeze(varcope_by_state(:,iContrast,iFreq,:)))
                             figure(5);imagesc(squeeze(first_level_results.stdcope(indind,:,iContrast,:))');title('stdcope');
                             figure(6);imagesc(squeeze(first_level_results.cope(indind,:,iContrast,:))');title('cope');
                             figure(7);imagesc(squeeze(first_level_results.cope(indind,:,iContrast,:)./first_level_results.stdcope(indind,:,iContrast,:))');title('tstat');
                             
                         end;
                    else
                        first_level_results.cope(indind,:,:,:)=first_level_results.cope_by_state(indind,:,:,:);
                        first_level_results.stdcope(indind,:,:,:)=sqrt(varcope_by_state);
                    end;
                    
                    first_level_results.stdcope_by_state(indind,:,:,:,:)=sqrt(varcope_by_state);
                                    
                    for iFreq = 1:size(first_level_results.cope_by_state,4)         
                         %% do bc on the copes
                        if(sum(first_level.bc)>0),
                            baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);

                            for c=1:length(contrast_list),
                                if(first_level.bc(c))
                                    first_level_results.cope(indind,:,c,iFreq)=first_level_results.cope(indind,:,c,iFreq)-mean(first_level_results.cope(indind,baseline_time_indices,c,iFreq),2);
                                end;
                            end;

                        end;
                    end;
                    
                else
                    first_level_results.cope(indind,:,1,:)=dat(:,:);
                    first_level_results.stdcope(indind,:,1,:)=1;                    
                    first_level_results.cope_by_state(indind,:,1,:)=dat(:,:);
                    first_level_results.stdcope_by_state(indind,:,1,:)=1;

                end;   
            else
                error('Outputting glm input data not implemented');
                if(~isfield(first_level_results,'glm_input_data')),
                    first_level_results.glm_input_times=glm_input_times;
                    first_level_results.glm_input_frequencies=out.tf_freqs;
                    first_level_results.glm_input_data=zeros([Nvoxels_out, size(dat)]);% num_voxels x  num_timepoints x num_contrasts x num_freqs                    
                end;                    
                    
                first_level_results.glm_input_data(indind,:,:,:)=dat;
            end;
                                                    
        end % if first_level.doGLM
    end % for indind = 1:Nvoxels_out
    
    if (first_level.doGLM && ~isempty(first_level.connectivity_seed_mni_coord)),
        first_level_results.cope_by_state(first_level_results.seed_index,:,:,:,:)=1;
        first_level_results.stdcope_by_state(first_level_results.seed_index,:,:,:,:)=1e16;
    end;
    
    ft_progress('close');   
    
    %% save results
    first_level_results.x=x;
    first_level_results.contrasts=first_level.contrast;
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.first_level=first_level;
    first_level_results.source_recon=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    try
        first_level_results.session_name=source_recon_results.session_name;
    catch,
        first_level_results.session_name=source_recon_results.subject_name;
    end;
    
    % save results to disk
    if first_level.save_data
        if ~isempty(first_level.trialwise_directory)
            trlWise_fname=fullfile(first_level.rawTF_directory, [first_level_results.session_name '_trlWise']);
        else
            first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
            trlWise_fname=fullfile(oat.source_recon.dirname, [first_level_results.session_name '_trlWise']);
        end
        disp(['Saving trialwise data in file ' oat.source_recon.dirname '/' first_level_results.fname]);
        save(trlWise_fname, 'trlwise', '-v7.3');
        clear trlwise
    end % if first_level.save_data
    
    if first_level.doGLM
        first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
        disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
        osl_save_oat_results( oat, first_level_results );
        results_fnames{sub}=first_level_results.fname;
    end % if first_level.doGLM
    
    if (strcmp(source_recon_results.recon_method,'none') | first_level.space_average),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
end;


